Modelling and Mitigating Secondary Crash Risk for Serial Tunnels on Freeway via Lighting-Related Microscopic Traffic Model with Inter-Lane Dependency

This paper models and mitigates the secondary crash (SC) risk for serial tunnels on the freeway which is incurred by traffic turbulence after primary crash (PC) occurrence and location-heterogeneous lighting conditions along serial tunnels. A traffic conflict approach is developed where SC risk is quantified using a surrogate safety measure based on the simulated vehicle trajectories after PC occurs from a lighting-related microscopic traffic model with inter-lane dependency. Numerical examples are presented to validate the model, illustrate SC risk pattern over time, and evaluate the countermeasures for SC, including adaptive tunnel lighting control (ATLC) and advanced speed and lane-changing guidance (ASLG) for connected vehicles (CVs). The results demonstrate that the tail of the stretching queue on the PC occurrence lane, the adjacent lane of the PC-incurred queue, and areas near tunnel portals are high-risk locations. In serial tunnels, creating a good lighting condition for drivers is more effective than advanced warnings in CVs to mitigate SC risk. Combined ATLC and ASLG is promising since ASLG informs CVs of an immediate response to traffic turbulence on the lane where PC occurs and ATLC alleviates SC risks on adjacent lanes via smoothing the lighting condition variations and reducing inter-lane dependency.


Introduction
In the last decade, the mileage of mountainous freeways has rapidly increased in the middle and western regions of China under the national freeway network improvement policy [1]. To improve the road alignment, tunnels are constructed and account for a significant proportion of mountainous freeways. Some of the tunnels are located close to each other and form a section of serial tunnels (successive tunnels, consecutive tunnels, tunnel group, equivalently) [2]. Tunnels are crash-prone in the daytime due to their narrowness with a semi-enclosed structure which leads to low illumination in the interior zone of the tunnel [3,4]. Drivers' behaviors are strongly influenced by lighting condition variations along serial tunnels, especially by the "black hole" and "white hole" effects near tunnel portals which frequently result in a crash risk [5]. To be specific, sharp lighting difference between the exterior and interior of the tunnel portal reduces the sight distance of the driver in the approaching vehicle, and then the driver experiences visual oscillation incurred by an abrupt lighting transition through the entrance and exit zones of the tunnel [6,7], as represented in Figure 1. In addition, once the crash occurs, a secondary crash (SC) is highly the sight distance of the driver in the approaching vehicle, and then the driver experiences visual oscillation incurred by an abrupt lighting transition through the entrance and exit zones of the tunnel [6,7], as represented in Figure 1. In addition, once the crash occurs, a secondary crash (SC) is highly likely to happen since drivers may not have immediate responses to the hazardous traffic turbulence [8]. Therefore, the frequent transitions between daylight and tunnel lighting in serial tunnels require drivers to visually adapt and react quickly, which increases the SC risk and challenges the safety condition of serial tunnels. SCs are crashes that occur within the spatial and temporal impact range of a primary crash (PC) [9][10][11]. As described above, SC risk in serial tunnels mainly results from traffic turbulence, which is caused by both PC occurrence and lighting variations along serial tunnels. PC-incurred traffic turbulence is a prevalently used contributor [12][13][14][15][16][17][18][19][20]. The impacts of traffic turbulence on SC risk are threefold. The first is the shockwave related to car-following behavior on the lane where PC occurs. The second is the vehicle type heterogeneity since trucks have a lower speed and reduced emergency brake ability in comparison with cars when they respond to the traffic turbulence [21][22][23]. The third is inter-lane dependency for multilane freeways, that vehicles' speeds are significantly reduced and a shockwave is formed by the queue on the adjacent lane, even though no vehicle in the queue attempts to change the lane [24][25][26]. However, the inter-lane dependency is hardly considered in SC risk modelling. In addition, to the best of our knowledge, there is no literature incorporating location-heterogeneous lighting conditions into the generation of traffic turbulence. Therefore, a lighting-related traffic model with inter-lane dependency is needed to reveal the mechanism of SC risk in serial tunnels.
Accordingly, the countermeasures for SC risk consist of two aspects: immediate informing of traffic turbulence and lighting variations smoothness along serial tunnels. Intelligent transportation systems are constantly being upgraded to make the countermeasures technically feasible. The emergence of connected vehicles (CVs) with advanced speed and lane-changing guidance (ASLG) is considered to be a promising solution for SC mitigation [27][28][29][30], even if its market penetration is at a low level up to date. Adaptive tunnel lighting control (ATLC) smooths the transition between daylight and tunnel lighting [31,32], such that the driver's visual oscillation with the vehicle's speed turbulence can be alleviated. However, whether ASLG and ATLC can effectively alleviate SC risk in serial tunnels needs to be further evaluated.
Hence, the aim of this study is to model the SC risk of serial tunnels on the freeway where the impacts of traffic turbulence after PC occurrence and lighting condition variations along serial tunnels on SC risk are described. The countermeasures for SC risk are evaluated, including ATLC for the present and ASLG for CVs in the future mobility.
The remainder of this paper is organized as follows. Section 2 provides a literature review on SC risk modelling and highlights our contribution. The model is formulated in Section 3. In Section 4, scenarios and simulations are conducted. Section 5 provides the SCs are crashes that occur within the spatial and temporal impact range of a primary crash (PC) [9][10][11]. As described above, SC risk in serial tunnels mainly results from traffic turbulence, which is caused by both PC occurrence and lighting variations along serial tunnels. PC-incurred traffic turbulence is a prevalently used contributor [12][13][14][15][16][17][18][19][20]. The impacts of traffic turbulence on SC risk are threefold. The first is the shockwave related to car-following behavior on the lane where PC occurs. The second is the vehicle type heterogeneity since trucks have a lower speed and reduced emergency brake ability in comparison with cars when they respond to the traffic turbulence [21][22][23]. The third is interlane dependency for multilane freeways, that vehicles' speeds are significantly reduced and a shockwave is formed by the queue on the adjacent lane, even though no vehicle in the queue attempts to change the lane [24][25][26]. However, the inter-lane dependency is hardly considered in SC risk modelling. In addition, to the best of our knowledge, there is no literature incorporating location-heterogeneous lighting conditions into the generation of traffic turbulence. Therefore, a lighting-related traffic model with inter-lane dependency is needed to reveal the mechanism of SC risk in serial tunnels.
Accordingly, the countermeasures for SC risk consist of two aspects: immediate informing of traffic turbulence and lighting variations smoothness along serial tunnels. Intelligent transportation systems are constantly being upgraded to make the countermeasures technically feasible. The emergence of connected vehicles (CVs) with advanced speed and lane-changing guidance (ASLG) is considered to be a promising solution for SC mitigation [27][28][29][30], even if its market penetration is at a low level up to date. Adaptive tunnel lighting control (ATLC) smooths the transition between daylight and tunnel lighting [31,32], such that the driver's visual oscillation with the vehicle's speed turbulence can be alleviated. However, whether ASLG and ATLC can effectively alleviate SC risk in serial tunnels needs to be further evaluated.
Hence, the aim of this study is to model the SC risk of serial tunnels on the freeway where the impacts of traffic turbulence after PC occurrence and lighting condition variations along serial tunnels on SC risk are described. The countermeasures for SC risk are evaluated, including ATLC for the present and ASLG for CVs in the future mobility.
The remainder of this paper is organized as follows. Section 2 provides a literature review on SC risk modelling and highlights our contribution. The model is formulated in Section 3. In Section 4, scenarios and simulations are conducted. Section 5 provides the results for the model validation, features of SC risks over time, and a countermeasure evaluation for SC risk. Conclusions of this paper and avenues for future research are presented in Section 6. Appendix A demonstrates the calibration of the lighting-related behavioral parameters in the model.

Literature Review
Compared with PCs which usually result from the driver's error, SCs are more deterministic and predictable because of their one-way causalities to PCs. The main contributing factors for SC risk include the type and duration of PC, time of day, road geometry, and traffic turbulence upstream of PC, where the last one is presented as the key factor [33].
Two mainstream approaches are developed to link the traffic turbulence with SC risk: crash-based statistical and traffic conflict approaches. In the former, SC risk is modelled as a direct measure, SC occurrence, described by a binary outcome (i.e., crash vs noncrash), according to the records in the crash database [33]. Traditional statistical models and machine learning techniques are developed to link the SC occurrence with traffic turbulence [12][13][14][15][16][17][18][19][20][34][35][36][37][38][39][40][41][42][43][44][45]. Traffic turbulence is quantified by macroscopic traffic flow variables along the spatial impact range of PC, such as average volume, speed, and occupancy, obtained from high-density probe vehicle data or traffic flow sensors with spacing of about 0.5 or 1.0 m on average.
However, this method is insufficient for serial tunnels on the Chinese freeway for two main reasons. The first is the unavailability of crash records and the lack of traffic flow sensors along the section. The second is that it cannot describe the traffic turbulence caused by lighting condition variations. The historical crashes are recorded in traffic police centers and are confidential to the public. The traffic-related data are collected from three types of cameras in serial tunnels installed by the freeway management center, as presented in Figure 2. The first is the camera at the overhead gantry a few kilometers upstream the first tunnel. It has the function of automatic vehicle identification (AVI) where the traffic flow variables are available. The second is the panoramic camera on the median of the freeway where a 1 km range can be covered. The third is the roadside camera with 150 m spacing on the tunnel wall which has a relatively narrow side view in the interior of tunnel. The latter two types of cameras merely have the function of incident detection, even though the views of them cover the entire serial tunnels. In other words, when they automatically capture the occurrence of PC, the traffic turbulence after PC occurs cannot be predicted since the real-time traffic flow variables along the serial tunnels are absent. results for the model validation, features of SC risks over time, and a countermeasure evaluation for SC risk. Conclusions of this paper and avenues for future research are presented in Section 6. Appendix A demonstrates the calibration of the lighting-related behavioral parameters in the model.

Literature Review
Compared with PCs which usually result from the driver's error, SCs are more deterministic and predictable because of their one-way causalities to PCs. The main contributing factors for SC risk include the type and duration of PC, time of day, road geometry, and traffic turbulence upstream of PC, where the last one is presented as the key factor [33].
Two mainstream approaches are developed to link the traffic turbulence with SC risk: crash-based statistical and traffic conflict approaches. In the former, SC risk is modelled as a direct measure, SC occurrence, described by a binary outcome (i.e., crash vs non-crash), according to the records in the crash database [33]. Traditional statistical models and machine learning techniques are developed to link the SC occurrence with traffic turbulence [12][13][14][15][16][17][18][19][20][34][35][36][37][38][39][40][41][42][43][44][45]. Traffic turbulence is quantified by macroscopic traffic flow variables along the spatial impact range of PC, such as average volume, speed, and occupancy, obtained from high-density probe vehicle data or traffic flow sensors with spacing of about 0.5 or 1.0 m on average.
However, this method is insufficient for serial tunnels on the Chinese freeway for two main reasons. The first is the unavailability of crash records and the lack of traffic flow sensors along the section. The second is that it cannot describe the traffic turbulence caused by lighting condition variations. The historical crashes are recorded in traffic police centers and are confidential to the public. The traffic-related data are collected from three types of cameras in serial tunnels installed by the freeway management center, as presented in Figure 2. The first is the camera at the overhead gantry a few kilometers upstream the first tunnel. It has the function of automatic vehicle identification (AVI) where the traffic flow variables are available. The second is the panoramic camera on the median of the freeway where a 1 km range can be covered. The third is the roadside camera with 150 m spacing on the tunnel wall which has a relatively narrow side view in the interior of tunnel. The latter two types of cameras merely have the function of incident detection, even though the views of them cover the entire serial tunnels. In other words, when they automatically capture the occurrence of PC, the traffic turbulence after PC occurs cannot be predicted since the real-time traffic flow variables along the serial tunnels are absent.  On the contrary, the traffic conflict approach is suitable for SC risk evaluation in serial tunnels on the Chinese freeway since it was initiated when the data for historical crashes were unavailable or traffic flow conditions were unknown. It can act as a surrogate for safety measures due to the causal relationship between traffic conflict severity and crash occurrence. The traffic conflict severity after PC occurrence is based on simulated vehicle trajectories from a microscopic traffic model and is evaluated by surrogate safety measures (SSMs) [46][47][48][49][50][51]. In microscopic traffic models, car-following (CF) and lane-changing (LC) maneuvers are described on a second or sub-second basis. For example, Li et al. [9] evaluated SC risk in adverse weather where historical crash data were scarce since adverse weather occurs infrequently. The simulated trajectories were obtained from a CF model with limited road surface adhesion and sight distance during adverse weather. The traffic conflict severity was quantified via time-to-collision (TTC), time exposed TTC, and time integrated TTC. Peng and Xu [11] evaluated the effectiveness of ASLG from arrays of roadside variable message signs on SC prevention where the evolutionary traffic condition under roadside ASLG after PC occurs can solely be predicted via simulation. The trajectories were obtained from the microscopic traffic simulator SUMO where CF and LC models. The traffic conflict severity was quantified via TTC and the deceleration rate to avoid the crash (DRAC). Son et al. [30] assessed in-vehicle ASLG in CV environment to prevent SCs. Since the driving behaviors are different in a CV environment, the trajectories were collected from driving simulation experiments and traffic conflict severity was quantified via TTC.
As described above, the key for traffic conflict-based SC analysis is the simulated vehicle trajectories after PC occurrence based on microscopic traffic models. Most CF and LC models are location-homogeneous where the behavioral parameters in the model remain constant during the course of driving [52]. However, they cannot accommodate vehicle maneuvers in serial tunnels where behavioral parameters are location-heterogeneous due to lighting condition variations along the section. In addition, the inter-lane dependency is hardly considered in traffic flow modelling to evaluate SC risk.
In terms of the abovementioned research gaps, the objective of this paper is to develop a traffic conflict approach where the SC risk in serial tunnels is quantified using SSM based on the vehicle trajectories after PC occurrence simulated from a lighting-related microscopic traffic model with inter-lane dependency (LMTID). The main contributions of this study are presented as follows: • A feasible solution is proposed for SC risk evaluation in serial tunnels of Chinese freeways where there is a lack of traffic flow sensors and an unavailability of crash records.

•
The behavioral parameters of CF and LC models are associated with the lighting conditions of different zones in serial tunnels. The link between SC risk in serial tunnels and lighting condition variation is modelled.

•
The inter-lane dependency is incorporated into the CF model to describe the impact of traffic turbulence from the PC-occurrence lane on the adjacent lane. The link between SC risk and traffic turbulence from PC is improved for the multilane freeway which conventionally considers the effects of vehicle type heterogeneity and CF-related shockwave on the lane where PC occurs.

Model Description
The key for traffic conflict-based crash risk analysis in serial tunnels is the simulated vehicle trajectories from a suitable microscopic traffic flow model. The lighting condition variations along the serial tunnels have an impact on drivers' desired speeds and perceived spacings. To be specific, (1) if the lighting difference (LD) between the exterior and interior of the tunnel portal is sharp, it is hard for the following vehicle to correctly perceive the spacing from preceding vehicles [7]; (2) if the lighting transition (LT) through the portal areas (i.e., entrance and exit zones) is abrupt, the driver has to experience the visual adaptation (VA) with slowing down the vehicle to accommodate the visual oscillation [53][54][55]; (3) the duration of VA reduces as driver travels through serial tunnels since he or she accommodates the consecutive LTs in the portal areas [56]; (4) due to the low illuminance in the interior zone, the driver keeps a relatively low desired speed with carefulness, even though he or she has lower discomfort after VA through the entrance zone [57,58].
Hence, the lighting-related desired speeds and perceived spacings in different locations should be presented by LMTID. A desired measures model is appropriate for the CF part of LMTID, which is one of the mainstream CF models [52]. Intelligent driver model (IDM) is one of the most popular time-continuous desired measures models using desired speed [59]. Due to its broad extrapolation, IDM is extended by following studies from various aspects [60][61][62][63][64][65]. Hence, it is feasible to incorporate the lighting-related desired speed and LD-related perceived spacing into IDM. In addition, the vehicle type heterogeneity and inter-lane dependency in CF can also be described by the IDM-formed model. The former has already been considered to reflect the varieties of speed and braking ability between cars and trucks' responses to the traffic turbulence [21][22][23]. The latter [24][25][26] has been incorporated into CF models, such as optimal velocity models [66,67], and the similar extension will be used in IDM.
The LC models are divided into two phases, LC decision-making (LCD) and LC process (LCP) modelling. In LCD, rule-based models appear to be the most popular [68][69][70] since they have deterministic rules for the probability of subject vehicle's LCD. The rule is mainly based on spacings between the subject and surrounding vehicles on multiple lanes. The LCD structure proposed by rule-based models is flexible and CF models can be easily integrated. However, in many of existing LC models, the second phase LCP modelling is not necessary and LC behavior is assumed to be executed instantaneously such that the lateral movement during LC and its effect on traffic are not clearly described. Actually, LCP has a significant impact on the following traffic on both the current and target lanes [71]. In addition, the travel speed during LCP is also responsive to the preceding vehicles on both the current and target lanes. Hence, CF models are used in LCP where the preceding vehicle for the subject vehicle is changeable [72].
It should be noted that behavioral decisions of the drivers might be determined by cultural differences [73,74]. Behavioral parameters, such as comfortable acceleration, CF spacing in a jam situation, the driver's reaction time, and the minimum necessary spacing are utilized from studies observed in China [23,75]. If LMTID is extended to other countries, the parameters need to be modified via field observations. However, the formulation of LMTID can remain unchanged.
In sum, an IDM-formed model is developed in the CF part of LMTID where lightingrelated desired speeds and perceived spacings, vehicle type heterogeneity, and inter-lane dependency are incorporated. The LC part of LMTID consists of a deterministic rule integrated with lighting-related perceived spacing in LCD, and the same CF model of LMTID in LCP.

Definitions, Notations and Assumptions
In a section of serial tunnels with I tunnels, give the ith tunnel with length D i . We use the road surface luminance as the approximation for driver visual luminance, as previous studies and current tunnel lighting standards did [76,77]. Let L(x) denote the road surface luminance (cd/m 2 ) with respect to the location x in the daytime where the beginning of the section is at x = 0, as demonstrated in Figure 2.
Serial tunnels are divided into four main zones with different characteristics of luminance in the daytime: exterior, entrance, interior, and exit zones, respectively, denoted as Z i o , Z i et , Z i in and Z i ex for the ith tunnel, respectively. The exterior zone is the open road outside where the exterior environmental luminance L 0 is constant. L 0 is the average luminance value focused in a conical aperture of a 20 • angle and looking forwards towards a centered point in the tunnel opening [77,78]. The entrance zone is the part of the tunnel where the lighting level is decreasing from the level at the end of the exterior zone to the level of the interior zone in the daytime. The interior zone is the inside of the tunnel when the luminance remains at a constantly low level, L i in . The exit zone is the part of the tunnel where the lighting level is increasing from the level at the end of the interior zone to the level of the exterior zone in the daytime. The entrance and exit zones are both referred to as portal areas. The lengths of exterior, entrance, interior and exit zones for the ith tunnel are denoted by D i o , D i et , D i in , and D i ex , respectively. The ending locations of exterior, entrance, interior, and exit zones for the ith tunnel are denoted by x i o , x i et , x i in , and x i ex , respectively. In current tunnel lighting standards [76,77], the zone divisions for serial tunnels are more sophisticated. The entrance zone is further divided into two subzones: the threshold and transition zones. The exterior zone is classified into three types: access, connecting, and parting zones. However, the simplified zone divisions in our study will not change the formulation of LMTID and can improve the readability of the paper.
The CF part of LMTID is developed via IDM. The lighting-related component for CF considers the desired speed of subject vehicle affected by LT in the portal area and low illuminance in the interior zone, and the perceived spacing of following vehicle affected by LD between the exterior and interior of the tunnel portal. As a result, the following definitions and assumptions are presented. Definition 1. The desired speed of the subject car at a location is the average speed at the same location when no surrounding car appears. Definition 2. LD for the subject vehicle is the ratio of luminance between the location of the preceding vehicle and location of the subject vehicle at the same time. If the spacing between the subject vehicle and the preceding vehicle is smaller than the perceived distance of the driver in the subject vehicle, the location and speed of the preceding vehicle can be correctly recognized by the subject vehicle. Otherwise, the preceding vehicle "disappears" in the vision of the driver in the subject vehicle, the spacing and speed difference perceived as maximum perceived distance and zero, respectively.

Assumption 2.
In the portal areas, whether the driver's VA happens is determined by LT. The VA duration and the desired speed reduction during VA are derived from LT, independently. If VA does not happen in the entrance zone, the desired speed linearly changes between the end of the exterior zone and the beginning of interior zone. If VA does not happen in the exit zone, the desired speed linearly changes between the end of interior zone and the beginning of exterior zone. Assumption 3. In serial tunnels, driver's VA duration is reduced by the number of times that VAs have been experienced.
Other assumptions in IDM are involved in the CF submodel of LMTID. Let x n (t), v n (t), and a n (t) denote the location, speed, and acceleration of the vehicle n at time t, respectively. For the vehicle n, the preceding vehicles on the same and adjacent lanes are denoted by the subscripts (n − 1) and (k − 1), respectively. The following vehicles on the same and adjacent lanes are denoted by the subscripts (n + 1) and k, respectively. ∆x 1 n (t) (∆x 2 n (t)) and s 1 n (t) (s 2 n (t)) denote the spacing from the preceding vehicles on the same (adjacent) lane to the vehicle n and the spacing from the front edge of the vehicle n to the rear end of the preceding vehicles on the same (adjacent) lane at time t, respectively. We have ∆x 1 where l n denotes the length of the vehicle n. The length varies according to vehicle type.
Let ∆v 1 n (t) (∆v 2 n (t)) denote the speed difference between the vehicle n and the preceding vehicle on the same (adjacent) lane, i.e., ∆v 1 ) denotes LD on the same (adjacent) lane at time t and we have L(x n (t)) ) according to Definition 2. L t (x n (t)) denotes LT at time t and we have L t (x n (t)) = L(x n (t)) L(x n (t−∆t)) where ∆t = 0.2 s according to Definition 3.

Lighting-Related Behavioral Parameters
Since the derivations of behavioral parameters in CF and LC models are similar, we introduce them in the CF scenario. When the driver accesses the tunnel portals, his or her perceived distance for the preceding vehicle is affected by LD. Let d 1 n (·) denote the perceived distance of the driver in the subject vehicle for the preceding vehicle on the same lane (m), derived from L 1 d (x n (t)). The LD-related d 1 n (·) has been investigated by previous studies via a gray target recognition experiment [7], as presented in Equation (1).
In Equation (1), the first and second equations show the perceived distances of drivers accessing entrance and exit portals, respectively. Letŝ 1 n (t) and ∆v 1 n (t) denote the perceived spacing and speed difference on the same lane, respectively. In terms of Assumption 1, we haveŝ 1 where s max is the maximum perceived spacing of drivers (m). When the driver travels through the portal areas, whether he or she experiences VA is affected by LT. Let T a (·) and b a (·) denote the VA duration (s) and the ratio of desired speed reduction during VA, respectively. Let k 1 and k 2 denote the thresholds for VA occurrences at entrance and exit zones, respectively. The LT-related T a (·) and b a (·) has been investigated by our field observation described in Appendix A, as presented in Equations (4) and (5): where k 1 = 0.025 and k 2 = 73, of which the derivations are described in Appendix A. When the driver traverses the serial tunnels, his or her VA duration is reduced by the number of times that VAs have been experienced according to Assumption 3. Let δ(ξ) denote the reduction coefficient on the VA duration where ξ denotes the number of times that VAs have been experienced. As the experimental results in Cui et al. [56] indicate, δ(ξ) is presented as the following.
Due to the constantly low illuminance in the interior zone, the driver keeps a constantly lower desired speed, even though he or she has lower discomfort after VA through the entrance zone. Let b in (·) denotes the ratio of the desired speed reduction in the interior zone with respect to the luminance L i in . The lighting-related b in (·) has been investigated by our field observation described in Appendix A, as presented in Equation (7).
In sum, the desired speeds V(·) for zones along serial tunnels are presented as follows. For exterior and interior zones, we have where v 0 n denotes the desired speed on the open road for the vehicle n. For portal areas, in terms of Assumption 2, if k 1 < L t (x n (t)) < k 2 , the VA will not happen and we have If L t (x n (t)) ≤ k 1 or L t (x n (t)) ≥ k 2 , the VA will happen and we have

CF Part of LMTID
Based on IDM, the dynamical equation of the vehicle n in the CF submodel of LMTID is given by where a max n and a comf n denote the maximum acceleration and comfortable deceleration, respectively. s jam n denotes the minimum spacing at the standstill situation. h n is the time gap. λ is the response coefficient of the relative speed between vehicles on the subject and adjacent lanes.
In the model, V(x n (t)) in Equation (11) presents the impact of location-heterogeneous lighting conditions along serial tunnels. The impact of traffic turbulence is reflected from three aspects. First, the model can simulate the shockwave after a given PC occurs in a recursive way. Second, the inter-lane dependency is described by the last term of Equation (12). Third, the different responses to traffic due to the vehicle type heterogeneity are presented by different behavioral parameters. For four types of vehicle-following patterns, desired speeds on the open road, minimum spacings at standstill situations, maximum accelerations, and comfortable accelerations are heterogeneous.
In terms of the dynamical equation, the key in the CF submodel of LMTID is the acceleration, a n (t), determined by five variables, v n (t), x n (t), s 1 n (t), ∆v 1 n (t), and ∆v 2 n (t). We denote acceleration as a n v n (t), x n (t),ŝ 1 n (t), ∆v 1 n (t), ∆v 2 n (t) for simplification.

LC Part of LMTID
The LC submodel of LMTID is divided into two phases, Lv et al.'s [75] incentive-based rule integrated with LD-related perceived spacing in LCD, and the same model of the CF submodel of LMTID in LCP. It should be noted that inter-lane dependency only exists in CF behavior while the queue forms on adjacent lane. Therefore, λ = 0 for LCP.
Incentive-based models are the variants of Gipps-type rule-based models which are governed by two LC rules, i.e., incentive and security criterion [79,80]. The idea behind the incentive-based models is intuitive and straightforward where there is only one "advantage" value, compared against a threshold value for final decision making. The subject and surrounding vehicles are demonstrated in Figure 3A. The rule in LCD is symmetrical and presented as follows.
Equation (16): where denotes the reaction time of the driver where = 1.45 s for drivers in cars and = 0.26 s for drivers in trucks [81]. and denote the friction coefficient and gravitational constant. Once the driver makes a decision to change lane, the LCP begins with a direction towards the target lane with an angle of , as demonstrated in Figure 3A. The formulation of is presented as follows: Vehicle n can change lane with a probability p 1 if Vehicle n can change lane with a probability p 2 if Vehicle n can change lane with a probability p 3 if Vehicle n does not change lane, otherwise. where n (t) and ∆v 2 n (t), respectively, of which the definitions are similar to Equations (2) and (3). s safe is the safe distance for lane changing where s safe = max s min , s b k . s min denotes the minimum necessary distance for the change. s b k denotes the braking distance of the following vehicle on the target lane (i.e., vehicle k in Figure 3A), as presented in Equation (16): where t r k denotes the reaction time of the driver where t r k = 1.45 s for drivers in cars and t r k = 0.26 s for drivers in trucks [81]. µ and g denote the friction coefficient and gravitational constant.
Once the driver makes a decision to change lane, the LCP begins with a direction towards the target lane with an angle of θ n , as demonstrated in Figure 3A. The formulation of θ n is presented as follows: where θ max and θ min represent the upper and lower bounds of the lane-changing angle, respectively. W denotes the lane width. Let y n (t) denote the lateral location of vehicle n at time t, where the middle position of the inner lane is the origin. After θ n is determined, the LCP is described in a recursive way, as shown in Table 1. During LCP, the inter-lane dependency is negligible such that λ = 0. In Algorithm 1, the scenario is taken as an example where vehicle n is changing from the inner to outer lane. Let lateral location of vehicle n in LCP y n (t) = 0, T denotes the time step.
Once the lane-changing of the subject vehicle is finished, it should be noted that vehiclefollowing patterns of the five vehicles (subject vehicle and four surrounding vehicles, see Figure 3) vary after LCP.

SSM-Based Crash Risk
The studies on SSM can date back to early 1970s [82] and considerable achievements have been made in this field since then. SSMs are typically estimated using the preceding and following vehicles' speeds and/or accelerations, and their spacing, which are derived from vehicle trajectory data based on recorded videos and/or microscopic traffic simulation results. SSMs fall into three categories: time-based SSMs, deceleration-based SSMs, and energy-based SSMs. The first two appear to be the most widely used. Time-based SSMs measure the time proximity to a crash occurrence, such as TTC [83]. Instead, decelerationbased SSMs focus on how vehicles' deceleration can prevent crash occurrence.
The most common deceleration-based SSM is DRAC, which was initially proposed by Cooper and Ferguson to determine the severity of a traffic conflict [84]. DRAC is the minimum deceleration rate required for a vehicle to avoid collision with a preceding vehicle. The formulation of DRAC is based on the assumption that the following vehicle takes evasive actions while the preceding vehicle retains its speed and direction. It was also extended to another form which is incorporated with driver's reaction times of the different types of following vehicles (i.e., car and truck) [81]. Compared to the time-based SSMs, the DRAC could reflect the evasive action of the following vehicle to stop in a timely manner or reach the corresponding speed of the leading vehicle and thereby avoid a crash [85]. Besides, the effectiveness of DRAC has been validated via comparison with the observed crash frequency [51,85].
The crash risk between two vehicles at the time R n (t) is defined as the probability that the instantaneous DRAC between two vehicles exceeds the maximum available deceleration rate (MADR), as Equation (19) indicates: where MADR n denotes MADR of vehicle n. MADR is associated with various factors, such as pavement condition and vehicle type. In this study, the MADR at the dry pavement condition is considered and is assumed to only depend on the vehicle type. MADR follows a truncated normal distribution and hence we have where Φ M  [81,86] are listed in Table 1.
DRAC is defined as the minimum deceleration rate of the following vehicle to stop in a timely manner behind the preceding vehicle as follows: Equation (27) gives the DRAC of vehicle n at time t in the CF scenario where the driver's reaction time t r n is considered. In the LC scenario, five vehicles (subject vehicle n and four surrounding vehicles in Figure 3B) are involved. Since DRAC describes the risk of following a vehicle colliding with the preceding vehicle, DRACs of vehicle n and its following vehicles on two lanes, vehicle n + 1 and vehicle k, are evaluated. Figure 3B describes the scenario where vehicle n is changing from the inner to the outer lane and we have DRAC n+1 (t) = DRAC (1) n+1 (t) y n ∈ −0.5W + w n 2 , y n+1 (t) + w n+1 2 + l n sin θ n + w n 2 cos θ n DRAC (2) n+1 (t) y n ∈ y n+1 (t) + w n+1 2 + l n sin θ n + w n 2 cos θ n , 1.5W − w n 2 where DRAC (1)

DRAC
(2) DRAC (1) n indicates the risk of vehicle n colliding with vehicle n − 1 while DRAC (2) n indicates the risk of vehicle n colliding with vehicle k. As vehicle n moves from inner to outer lane, the crash target of vehicle n transfers from vehicle n − 1 to vehicle k. If the lateral spacing between vehicle n − 1 and vehicle k is wide, there will be no crash risk for vehicle n at a certain time, as Equation (28) presents. If the lateral spacing between vehicle n − 1 and vehicle k is narrow, there will be a risk for vehicle n to collide with both vehicle n − 1 and vehicle k at a certain time, as Equation (29) presents. The core of the derivation of DRAC is similar while vehicle n is changing from the outer to inner lane and is not explicated in this paper.

Scenarios and Simulations
In this section, real-world-based numerical examples are presented from a section of serial tunnels on the two-lane freeway G65, Chongqing, China, as demonstrated in 2000 m, 2 = 2600 m, and 3 = 1200 m, respectively. The distances between two tunnels are 2 = 710 m and 3 = 780 m, respectively. The road is designed with solid lane marking stretching from 150 m ahead of the entrance portal of the 1st tunnel to 100 m behind of the exit portal of the 3rd tunnel, where LCs are prohibited for drivers to make. In other parts of serial tunnels, LCs are discretionary with dash lane marking. The widths for both lanes are 3.75 m. The noon or afternoon is the time for the basic scenario where exterior environmental luminance 0 = 6000 cd m 2 ⁄ . In the tunnel lighting design standard of China [77], a stepwise curve method is used for zone luminance. For the th tunnel the entrance zone is further divided into two portions of the threshold zone ℤ ℎ1 and ℤ ℎ2 , and three portions of the transition zone ℤ 1 , ℤ 2 , and ℤ 3 . The exit zone is further divided into two portions ℤ 1 and ℤ 2 . The ending positions of these portions are ℎ1 , ℎ2 , 1 , 2 , 3 , 1 , and 2 , respectively, where 3 = and 2 = . The lengths of these portions are

Scenarios
There are six possible locations of PC occurrences investigated, consisting of middle positions in all three tunnels on both inner and outer lanes, as shown in Figure 5.

Scenarios
There are six possible locations of PC occurrences investigated, consisting of middle positions in all three tunnels on both inner and outer lanes, as shown in Figure 5.
2000 m, 2 = 2600 m, and 3 = 1200 m, respectively. The distances between two tunnels are 2 = 710 m and 3 = 780 m, respectively. The road is designed with solid lane marking stretching from 150 m ahead of the entrance portal of the 1st tunnel to 100 m behind of the exit portal of the 3rd tunnel, where LCs are prohibited for drivers to make. In other parts of serial tunnels, LCs are discretionary with dash lane marking. The widths for both lanes are 3.75 m. The noon or afternoon is the time for the basic scenario where exterior environmental luminance 0 = 6000 cd m 2 ⁄ . In the tunnel lighting design standard of China [77], a stepwise curve method is used for zone luminance. For the th tunnel the entrance zone is further divided into two portions of the threshold zone ℤ ℎ1 and ℤ ℎ2 , and three portions of the transition zone ℤ 1 , ℤ 2 , and ℤ 3 . The exit zone is further divided into two portions ℤ 1 and ℤ 2 . The ending positions of these portions are ℎ1 , ℎ2 , 1 , 2 , 3 , 1 , and 2 , respectively, where 3 = and 2 = . The lengths of these portions are ℎ1 , ℎ2 , 1 , 2 , 3 , 1 , and 2 , respectively. The location at 1000 m ahead of the entrance portal of the 1st tunnel is denoted as = 0. Hence, 1 = 1 = 1000 m.

Scenarios
There are six possible locations of PC occurrences investigated, consisting of middle positions in all three tunnels on both inner and outer lanes, as shown in Figure 5.  An amount of 3000 cd/m 2 of the exterior environmental luminance L 0 is always the fundamental parameter of tunnel lighting design in China for all portions [77]. The two basic luminance levels are L i th1 = 0.025 × 3000 = 75 cd/m 2 and L i in = 1.0 cd/m 2 . The "black hole" and "white hole" effects are very likely to happen for drivers in the basic scenario where L 0 = 6000 cd/m 2 , as presented in Scenario 0.

Simulation Time Step Comparison
We determined the appropriate time step in the simulation for the sake of efficiency based on Scenario 0. The alternative time steps are 0.1 and 1.0 s. To make a fast decision on the time step, we reduce the scale of the car platoon and the duration of PC. In the simulation, 200 vehicles with 20% trucks run from the beginning of the section. The average values of initial spacings and desired speeds are 100 m and 80 km/h, respectively. The PC occurs on the outer lane of the middle positions in the 2nd tunnel (i.e., Location 2 in Figure 5). The duration for PC is set as 10 min.

Model Validation
The LMTID-based simulation is validated via the comparison with the experimental observation of real traffic in studied serial tunnels. In the simulation, the location and duration of PC occurrence are set in terms of the video records. In addition, the video records when no PC occurs are used for validation as well.
The CF model is validated in the area where the road is designed with solid lane marking and LCs are prohibited for drivers to make, as presented in Figure 5. Since the sequence of cars and the result of the CF model are deterministic, a platoon-based comparison between real traffic and a simulation is feasible, such as the distribution of travel times and RMSE of the trajectory's speed profile [87] in car platoon. The real trajectories of cars are constructed via the array of roadside cameras with 150 m spacing on the tunnel wall (third type of camera in Figure 2). Car speeds and times crossing the camera view (regarded as a virtual cross-sectional detector) are automatically captured by each camera. The cars' number plates and types (i.e., car or truck) are manually recorded in the video replay of each camera due to the narrow side of the camera and low illuminance in the interior zone of the tunnel. Based on the sequences of matched cars' number plates, types, times, and speeds crossing the virtual cross-sectional detectors in all cameras, trajectories of car platoons can be formed. In the simulation for the CF model, the luminance levels along serial tunnels should have been measured in field. However, we did not grant the permission to measure the road surface luminance in the studied serial tunnels. Therefore, we use the designed tunnel lighting scheme (i.e., Scenario 0) to substitute for the real-world road surface luminance. In addition, the initial state of the simulation coincides with that in the beginning of the video records, including the locations, speeds, and types of all vehicles.
The LC model is validated in the area at two ends of the section of serial tunnels where LCs are discretionary with dash lane marking and covered by a panoramic camera (second type of camera in Figure 2), as presented in Figure 5. Due to the stochasticity of LCs, platoon-based validation is insufficient. Instead, we use macroscopic traffic variables of real traffic and a simulation for comparison, such as a density-flow diagram and LC frequency [75]. The experimental observations are recorded by panoramic camera. The flow and density are derived from the vehicle number of each lane for every second. In the simulation for the LC model, the initial state of the simulation coincides with that in the beginning of the video records, including the locations, speeds, and types of all vehicles.
Besides, the lighting-related behavioral parameters (T a , b a , k 1 , k 2 , and b in ) are calibrated via naturalistic driving experiments for drivers in other different tunnels, which are explicated in Appendix A.

SC Risk under Current Situation and Countermeasures
Firstly, we investigate the features of SC risks over time for different locations of PC occurrence in serial tunnels. Therefore, PC occurs in parallel at all six locations of Figure 5 in Scenario 0.
Then, the effectiveness of countermeasures is investigated. For the scenario under ATLC, the luminance levels in all portions adapt to the exterior environmental luminance 6000 cd/m 2 , L i th1 = 0.025 × 6000 = 150 cd/m 2 , L i in = 2.5 cd/m 2 [77], as presented in Scenario 1. For the scenario under ASLG, CVs will start to decelerate 1 km ahead of the tail of the queue if they decide to select the lane where PC occurs. The desired speed V(·) is assumed to be reduced by 40%. Otherwise, CVs will change to the adjacent lane 1 km ahead of the tail of the queue or 1 km ahead of the beginning of the solid lane marking. The proportions of CVs are 30%, 60%, and 90% for Scenarios 2, 3, and 4, respectively. Finally, the countermeasure combining ATLC and ASLG is proposed with 30% CVs, as Scenario 5 presents. In these scenarios, the locations of PC occurrences are the same, the inner lane of middle positions in 2nd tunnel (i.e., Location 2 in Figure 5).  In the simulations, 2500 vehicles with 20% trucks run from the beginning of the section. The average values of initial spacings and desired speeds are 100 m and 80 km/h, respectively. The duration for PC is set as 60 min.

Simulations
In the simulation, the vehicle lengths for cars and trucks are 6 m and 12 m, respectively. The vehicle widths for cars and trucks are 1.8 m and 2.5 m, respectively. For the CF submodel and LCP in LMTID, The response coefficient of the relative speed between vehicles on the subject and adjacent lanes λ = 0.6. The acceleration exponent β needs to be calibrated. For LCD in LMTID, the reaction times of the drivers in cars and trucks are 1.45 s and 0.26 s, respectively. The upper and lower bounds of the LC angle are θ min = 5 • and θ max = 20 • , respectively. The friction coefficient µ = 0.8 and the gravitational constant g = 9.8 m/s 2 . The LC probabilities [p 1 , p 2 , p 3 ] are the parameters that need to be calibrated. The minimum necessary spacing for LC s min = 5 m. For SC risk evaluation, the distribution parameters of MADR for cars and trucks are listed in Algorithm 1. The simulation is coded in MATLAB and run on a workstation.
For model validation, the simulation program is initialized according to the beginning frame of the video records, including the locations, speeds, and types of all vehicles. For the time step determination and SC risk evaluation, the simulation program is initialized by distributing vehicles randomly on the two lanes via the distribution of the spacings between two consecutive vehicles. The spacings and speeds yield to normal distribution. Whether a vehicle is a car or a truck, connected or unconnected vehicles are generated via Bernoulli distributions in terms of proportions, where no truck is generated on the inner lane initially.
In the program, we set an attribute set for each vehicle which includes vehicle No., lane No., serial No., vehicle angle, position, velocity, acceleration, and spacing. The vehicle's serial No. remains unchangeable while other attributes will be updated during the course of simulation. In each time step, a CF calculation (i.e., Equations (11) and (12)) is accomplished to determine the vehicle movement on the subject lane. Then a judgment for LCD is made, whether a vehicle changes lanes or not (i.e., Equations (13) and (15)). If the LC rule is satisfied, the vehicle will start LCP and LCP is operated according to Algorithm 1. At the end of LCP, the vehicle has moved to the target lane and then it will move on adjacent lane in the next time step. After each time step, all vehicles' attributes, i.e., lane, angle, position, velocity, acceleration, and spacing, are updated.
It should be noted that the initial state of the simulation can be set according to the real-time traffic flow data from the AVI camera (first type of camera in Figure 2) if the program is plugged into the traffic surveillance system in the freeway management center.

Simulation Time Step Determination
Since the vehicle trajectories are stochastic from the simulation, we cannot directly compare the trajectories between two time steps. Instead, we use the mean speed in a spatiotemporal interval as the comparison metric. The temporal and spatial intervals are 2 min and 200 m, respectively. After 100 parallel simulations, the medians of metrics in the two time steps (i.e., 0.1 and 1.0 s) and their difference are presented as contours in Figures 6 and 7. hicle's serial No. remains unchangeable while other attributes will be updated during the course of simulation. In each time step, a CF calculation (i.e., Equations (11) and (12)) is accomplished to determine the vehicle movement on the subject lane. Then a judgment for LCD is made, whether a vehicle changes lanes or not (i.e., Equations (13) and (15)). If the LC rule is satisfied, the vehicle will start LCP and LCP is operated according to Algorithm 1. At the end of LCP, the vehicle has moved to the target lane and then it will move on adjacent lane in the next time step. After each time step, all vehicles' attributes, i.e., lane, angle, position, velocity, acceleration, and spacing, are updated.
It should be noted that the initial state of the simulation can be set according to the real-time traffic flow data from the AVI camera (first type of camera in Figure 2) if the program is plugged into the traffic surveillance system in the freeway management center.

Simulation Time Step Determination
Since the vehicle trajectories are stochastic from the simulation, we cannot directly compare the trajectories between two time steps. Instead, we use the mean speed in a spatiotemporal interval as the comparison metric. The temporal and spatial intervals are 2 min and 200 m, respectively. After 100 parallel simulations, the medians of metrics in the two time steps (i.e., 0.1 and 1.0 s) and their difference are presented as contours in Figures 6 and 7.  As Figures 6 and 7 indicate, the predicted traffic flow evolutions after PC occurs with different time steps are similar. However, the average running times in 100 parallel simulations with 0.1-and 1.0-min time steps are 5.4 and 0.6 min, respectively. As a result of a lower computational capacity requirement and a higher efficiency, 1.0 s is determined as the time step for the simulation in all following scenarios. In addition, the relatively faster speed of the program with a 1.0 s time step makes our approach more feasible to plug into the traffic surveillance system in the freeway management center. As Figures 6 and 7 indicate, the predicted traffic flow evolutions after PC occurs with different time steps are similar. However, the average running times in 100 parallel simulations with 0.1-and 1.0-min time steps are 5.4 and 0.6 min, respectively. As a result of a lower computational capacity requirement and a higher efficiency, 1.0 s is determined as the time step for the simulation in all following scenarios. In addition, the relatively faster speed of the program with a 1.0 s time step makes our approach more feasible to plug into the traffic surveillance system in the freeway management center.

Model Validation
For the CF part of LMTID, there are many parameters that need to be calibrated. Due to limited observation data in our study, the most sensitive parameter is selected to calibrate, and the values of others come from the results of previous studies observed in China. The acceleration exponent β is often presented as the most sensitive parameter to affect the results [59][60][61] and will be validated in this paper. The alternative values for β are 2 and 4. The values of h n , s jam n , v 0 n , a max n , and a comf n for four types of vehicle-following patterns are listed in Table 2 [23]. A 30 min video record after a PC occurs is available for us to validate the CF model. The PC was a slight rear-end crash between two cars, and hence the duration of PC lasted for about 18.4 min. The PC is located on the outer lane and 784 m to the entrance of the first tunnel. We compare the distributions of travel times and the RMSE of the trajectory's speed profile in car platoons on both lanes. The latter measure describes the magnitude of an individual vehicle's speed oscillation. The calculation of RMSE of the trajectory's speed profile can be found in [87]. The comparisons among observations and simulations with different values of acceleration exponent β are presented in Figure 8.   The results show that = 4 is appropriate for the CF model and the CF model roughly corresponds with the experimental observation result. The error may result from the narrowness effect of the tunnel wall and visual fatigue due to low illuminance [88] when drivers travel through serial tunnels. These factors affecting driving behavior need to be considered in LMTID in future research. In addition, distributions of travel times and the RMSE of the trajectory's speed profile in car platoons are bimodal where the second peak presents the vehicles congested upstream of the PC.
For the LC part of LMTID, LC probability is the most representative parameter that needs to be calibrated. The same video record with the CF model validation is used, with other video records when no crash occurs supplemented. We compare the flow-density diagram and LC frequency on both lanes. Since the results from the LC model are stochastic, we use the distribution of 100 parallel simulation results or the median of the distribution. The comparisons among observation and simulations with [ ] The results show that β = 4 is appropriate for the CF model and the CF model roughly corresponds with the experimental observation result. The error may result from the narrowness effect of the tunnel wall and visual fatigue due to low illuminance [88] when drivers travel through serial tunnels. These factors affecting driving behavior need to be considered in LMTID in future research. In addition, distributions of travel times and the RMSE of the trajectory's speed profile in car platoons are bimodal where the second peak presents the vehicles congested upstream of the PC.
For the LC part of LMTID, LC probability is the most representative parameter that needs to be calibrated. The same video record with the CF model validation is used, with other video records when no crash occurs supplemented. We compare the flowdensity diagram and LC frequency on both lanes. Since the results from the LC model are stochastic, we use the distribution of 100 parallel simulation results or the median of the distribution. The comparisons among observation and simulations with different values of LC probabilities [p 1 , p 2 , p 3 ] are presented in Figure 9. In the flow-density comparison, we compare the median of 100 parallel simulation results with the observation. [ 1 , 2 , 3 ] = [0.3,0.1,0.05] is consistent with the experimental observation. The results also imply that some of the drivers seem somewhat aggressive when changing lanes, which coincides with the driving characteristics of Chinese drivers in the studies of driving cultures [73,74].
In the LC frequency comparison, we compare the distribution of 100 parallel simulation results with the observation. Figure 9A demonstrates the results of the simulation agree well with the observation. As traffic density increases, LC frequency first increases and then decreases. This is because that vehicles change lanes for better driving conditions with the increase of traffic density and they find it difficult to change lanes with a small spacing to the surrounding vehicles on both lanes at a high density condition.
It should be noted that, to obtain more accurate results, more video records after PC occurs are needed to calibrate the LMTID. When the video record covering both PC and SC becomes available, DRAC-based SC risk in our study will be verified.

Features of SC Risks over Time for Different PC Locations in Serial Tunnels
For Scenario 0, we investigate the features of SC risks over time for different PC occurrence locations in serial tunnels. The high-risk points (HRPs) that are spatiotemporally distributed are aggregated for a better illustration of crash risk evolutions. A HRP is defined as ( , ( )) when ( ) ≥ 0.8 for vehicle . The temporal and spatial intervals for aggregation are 2 min and 200 m, respectively. Then, the density of HRPs can be obtained via the ratio of the aggregated number of HRPs to the area of the spatiotemporal interval. Traffic flow and crash risk evolutions for two lanes are presented in The results also imply that some of the drivers seem somewhat aggressive when changing lanes, which coincides with the driving characteristics of Chinese drivers in the studies of driving cultures [73,74].
In the LC frequency comparison, we compare the distribution of 100 parallel simulation results with the observation. Figure 9A demonstrates the results of the simulation agree well with the observation. As traffic density increases, LC frequency first increases and then decreases. This is because that vehicles change lanes for better driving conditions with the increase of traffic density and they find it difficult to change lanes with a small spacing to the surrounding vehicles on both lanes at a high density condition.
It should be noted that, to obtain more accurate results, more video records after PC occurs are needed to calibrate the LMTID. When the video record covering both PC and SC becomes available, DRAC-based SC risk in our study will be verified.

Features of SC Risks over Time for Different PC Locations in Serial Tunnels
For Scenario 0, we investigate the features of SC risks over time for different PC occurrence locations in serial tunnels. The high-risk points (HRPs) that are spatiotemporally distributed are aggregated for a better illustration of crash risk evolutions. A HRP is defined as (t, x n (t)) when R n (t) ≥ 0.8 for vehicle n. The temporal and spatial intervals for aggregation are 2 min and 200 m, respectively. Then, the density of HRPs can be obtained via the ratio of the aggregated number of HRPs to the area of the spatiotemporal interval. Traffic flow and crash risk evolutions for two lanes are presented in Figures 10-15 when PC occurs on both lanes at six locations, respectively. occurrence locations in serial tunnels. The high-risk points (HRPs) that are spatiotemporally distributed are aggregated for a better illustration of crash risk evolutions. A HRP is defined as ( , ( )) when ( ) ≥ 0.8 for vehicle . The temporal and spatial intervals for aggregation are 2 min and 200 m, respectively. Then, the density of HRPs can be obtained via the ratio of the aggregated number of HRPs to the area of the spatiotemporal interval. Traffic flow and crash risk evolutions for two lanes are presented in     As Figures 10-15 indicate, we can observe the followings.
• The stretching tail of the queue incurred by PC occurrence is always the high risky location over time, since the arriving vehicle has to brake to a halt before colliding with the tail of the queue.

•
For the lane adjacent to the PC occurrence lane, it has even higher numbers of high-risk points due to the inter-lane dependency and low illumination inside the tunnel. The spatial range of HRPs stretches as the PC-incurred queue on the adjacent lane extends. As a result of the "black hole" and "white hole" effects, the density of HRPs is significantly high near tunnel portals in serial tunnels.

•
If PC occurs at the rear location of serial tunnels, the density of SC risk on the PC occurrence lane is relatively low since drivers maintain a slower speed. However, the density of SC risk on the adjacent lane is higher with the stretching PC-incurred As Figures 10-15 indicate, we can observe the followings.
• The stretching tail of the queue incurred by PC occurrence is always the high risky location over time, since the arriving vehicle has to brake to a halt before colliding with the tail of the queue.

•
For the lane adjacent to the PC occurrence lane, it has even higher numbers of high-risk points due to the inter-lane dependency and low illumination inside the tunnel. The spatial range of HRPs stretches as the PC-incurred queue on the adjacent lane extends. As a result of the "black hole" and "white hole" effects, the density of HRPs is significantly high near tunnel portals in serial tunnels.

•
If PC occurs at the rear location of serial tunnels, the density of SC risk on the PC occurrence lane is relatively low since drivers maintain a slower speed. However, the density of SC risk on the adjacent lane is higher with the stretching PC-incurred As Figures 10-15 indicate, we can observe the followings.
• The stretching tail of the queue incurred by PC occurrence is always the high risky location over time, since the arriving vehicle has to brake to a halt before colliding with the tail of the queue.

•
For the lane adjacent to the PC occurrence lane, it has even higher numbers of highrisk points due to the inter-lane dependency and low illumination inside the tunnel. The spatial range of HRPs stretches as the PC-incurred queue on the adjacent lane extends. As a result of the "black hole" and "white hole" effects, the density of HRPs is significantly high near tunnel portals in serial tunnels.

•
If PC occurs at the rear location of serial tunnels, the density of SC risk on the PC occurrence lane is relatively low since drivers maintain a slower speed. However, the density of SC risk on the adjacent lane is higher with the stretching PC-incurred queue along serial tunnels, due to the inter-lane dependency and low illumination. • Periodic oscillations of the traffic flow pattern are presented near tunnel portals where VAs are needed for drivers, demonstrated by the high density traffic flow departing from the PC occurrence location after PC removal.

Comparison of Countermeasures for SC Risk
Traffic flow and crash risk evolutions are compared with different countermeasures for SC risk from Scenarios 1 to 5 when PC occurs on the inner lane of Location 2, as presented in Figures 16-20. Meanwhile, the effectiveness of countermeasures is compared with the scenario under no countermeasure, i.e., Figure 12.    As Figures 12 and 16-20 indicate, we can observe the followings.
• In serial tunnels, creating good lighting conditions for drivers is more effective than advanced warnings for a part of drivers in CVs to mitigate the SC risk after PC occurrence.

•
Since ATLC can successfully eliminate the "black hole" and "white hole" effects in terms of the daylight theoretically, the SC risk at the tail of the PC-incurred queue is alleviated when tail reaches the portal areas.  As Figures 12 and 16-20 indicate, we can observe the followings.
• In serial tunnels, creating good lighting conditions for drivers is more effective than advanced warnings for a part of drivers in CVs to mitigate the SC risk after PC occurrence.

•
Since ATLC can successfully eliminate the "black hole" and "white hole" effects in terms of the daylight theoretically, the SC risk at the tail of the PC-incurred queue is alleviated when tail reaches the portal areas.   A countermeasure combining ATLC with ASLG is promising since it can prevent the "black hole" and "white hole" effects near tunnel portals, inter-lane dependency, and low illumination inside the tunnel, and immediately inform drivers in CVs of the hazardous traffic flow evolution from PC.

Conclusions
To model and mitigate the secondary crash (SC) risk for serial tunnels on the freeway, this paper developed a traffic conflict approach where traffic turbulence is associated with primary crash (PC) occurrence and lighting condition variations along serial tunnels. In the approach, SC risk was quantified using surrogate safety measure (SSM) based on simulated vehicle trajectories after PC occurs via a lighting-related microscopic traffic model with inter-lane dependency (LMTID). CF and LC behaviors in LMTID were associated with location-heterogeneous lighting conditions along serial tunnels and traffic turbulence, which consists of shockwave after PC occurs, vehicle type heterogeneity, and inter-lane dependency. The deceleration rate to avoid the crash (DRAC) was utilized as a SSM to evaluate the SC risk. The numerical examples of a section of serial tunnels on a two-lane freeway were presented to validate the model, illustrate the SC risks over time for different PC locations and countermeasures, including adaptive tunnel lighting control (ATLC) for the present and advanced speed and lane-changing guidance (ASLG) for connected vehicles (CVs) in the future mobility.
The result demonstrates that the tail of the stretching queue on PC occurrence lane, the stretching section on adjacent lane of the PC-incurred queue, and areas near tunnel portals are the high-risk locations, due to the lack of immediate response to the hazardous traffic flow evolution from PC, inter-lane dependency, and low illumination inside the tunnel, and the "black hole" and "white hole" effects, respectively. If PC occurs at the rear location of serial tunnels, the density of high SC risk on the PC occurrence lane is relatively low since drivers maintain a slower speed when traversing the serial tunnels. Periodic oscillations of the traffic flow pattern are presented near tunnel portals when traffic flow departs from the PC occurrence location after PC removal.
In serial tunnels, creating a good lighting condition for drivers is more effective than advanced warnings in CVs to mitigate the SC risk after PC occurrence. To be specific, ATLC can alleviate the SC risk at the tail of a PC-incurred queue when the tail reaches the portal areas through eliminating the "black hole" and "white hole" effects. ATLC can considerably reduce the impact of inter-lane dependency and low illumination on crash risk along the adjacent lane of PC occurrence. ASLG for CVs manages to better alleviate the SC risk at the tail of the PC-incurred queue via immediately informing drivers in CVs of the hazardous traffic flow evolution from PC. However, a few new SC risky locations appear which result from the advance deceleration of CVs whose drivers decide not to change to the adjacent lane. Additionally, ASLG fails to limit the crash risk along the adjacent lane of PC occurrence since it cannot eliminate the inter-lane dependency after it guides a part of CVs to the adjacent lane. Hence, countermeasure combining ATLC with ASLG is promising.
Further research directions can be identified from the limitations of this paper. More factors affecting driving behavior need to be considered in LMTID, such as the narrowness effect of tunnel wall and visual fatigue travelling through serial tunnels. More video records after PC occurs are needed to more sophisticatedly calibrate and verify the LMTID and SC risk.  tunnel with its particular lighting condition. Points are presented only if VA occurs (i.e., VA duration is bigger than 0 s).
for all 50 drivers. The measurements of position, illuminance, car's travel speed, and pupil area of drivers were tuned to ensure the consistency of the timeline of measurements.

A.1. VA Duration and Desired Speed Reduction in Portal Areas
In the analysis, the average value for VA durations is used if the values of some LT are similar in different days or tunnels. The relationship between and is demonstrated in Figure A1. In Figure A1, each point (i.e., sample) is derived from a driver traversing a tunnel with its particular lighting condition. Points are presented only if VA occurs (i.e., VA duration is bigger than 0 s). As Figure A1 indicates, the logarithm function is considered to be appropriate for quantifying the impact of LT on VA duration, as presented in Equation (A1).
T a (L t (x n (t))) = −2.976 ln L t (x n (t)) − 11.188(s) L t (x n (t)) ≤ k 1 , R 2 = 0.78 1.275 ln L t (x n (t)) − 5.470(s) L t (x n (t)) ≥ k 2 , R 2 = 0.77 , where R denotes the coefficient of determination reflecting the goodness of curve fit. In addition, the threshold for whether the VA occurs is determined as the intercept on the vertical coordinate of the logarithm function. The thresholds for dark and bright adaptations are k 1 = 0.025 and k 2 = 73, respectively. In the desired speed reduction analysis, for each period (morning or afternoon) of a particular tunnel, the speed data of all 25 drivers is considered as a sample. The data of drivers are removed where VA does not occur. The relationship between L t and b a is demonstrated as Figure A2.
Linear function is appropriate to describe the impact of LT on the ratio of desired speed reduction, as presented in Equation (A2). As Figure A1 indicates, the logarithm function is considered to be appropriate for quantifying the impact of LT on VA duration, as presented in Equation (A1).
where denotes the coefficient of determination reflecting the goodness of curve fit. In addition, the threshold for whether the VA occurs is determined as the intercept on the vertical coordinate of the logarithm function. The thresholds for dark and bright adaptations are 1 = 0.025 and 2 = 73, respectively.
In the desired speed reduction analysis, for each period (morning or afternoon) of a particular tunnel, the speed data of all 25 drivers is considered as a sample. The data of drivers are removed where VA does not occur. The relationship between and is demonstrated as Figure A2. Linear function is appropriate to describe the impact of LT on the ratio of desired speed reduction, as presented in Equation (A2).

A.2. Desired Speed Reduction in Interior Zone
The desired speed of a particular lighting condition in the interior zone is determined as the average value of collected speeds at all locations where no preceding vehicle appears in the visions of drivers. For each driver, the average value for desired speed b a (L t (x n (t))) = 41.841L t (x n (t)) − 0.059 L t (x n (t)) ≤ k 1 , R 2 = 0.93 −0.001L t (x n (t)) + 1.100 L t (x n (t)) ≥ k 2 , R 2 = 0.96 (A2)